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In this article, a unified approach to obtain symplectic integrators on T*G from 
Lie group integrators on a Lie group G is presented. The approach is worked out 
in detail for symplectic integrators based on Runge-Kutta-Munthe-Kaas methods 
and Crouch-Grossman methods. In both cases, we show that it is possible to 
obtain symplectic integrators of arbitrarily high order by this approach. 

1 Introduction 

1.1 Motivation and background 

In general, an ordinary differential equation (ODE) can be described by a vector field on a 
smooth manifold where solutions of the ODE are integral curves of the vector field. Numerical 
approximation of solutions of ODEs is an old field of study, and a plethora of methods for 
obtaining numerical solutions exist. However, most of these methods assume that the manifold 
is Euclidean space. If the manifold is not Euclidean space, it is possible to embed the manifold 
in Euclidean space, and extend the vector field on the manifold to a vector field in Euclidean 
space such that the integral curves are ensured to remain in the image of the embedding. A 
standard numerical algorithm will in general result in discrete points which do not lie in the 
image of the embedding. An improvement of this approach is to use projection methods to 
obtain solutions on the manifold. These approaches, though simple, suffer from the problem 
that the numerical solutions depend on the particular choice of embedding, and on the 
particular extension of the vector field. 

One aspect of geometric numerical integration is to exploit structure on the manifold to 
define numerical methods that are intrinsic to the manifold (i.e. do not depend on a particular 
embedding). This structure can for instance be that of a Lie group acting on the manifold. The 
action of a Lie group G on a manifold M is a smooth mapping W:G x M -> M which respects 
the group structure on G. If the action is transitive, then the derivative with respect to the first 
component of f at the group identity e is a surjective vector bundle morphism qx M— TM. 
Any vector field X on M can then be lifted (possibly not in a unique manner) to a section of 
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the vector bundle q x M. The combination of this lifting and standard charts g—>G, form the 
basis of several classes of Lie group methods. Among them are the Crouch-Grossman (CG) 
methods |4| and the Runge-Kutta-Munthe-Kaas (RKMK) methods (12). For a more detailed 
discussion of Lie group methods, we refer to the survey article by Iserles et al. |8| and the 
references therein. 

Another aspect of geometric numerical integration is symplecticity of numerical integrators. 
Many important problems from physics can be formulated as Hamiltonian ODEs on cotangent 
bundles over manifolds. The flow maps of these ODEs are symplectic, that is, they preserve 
the canonical two -form on the cotangent bundle. For symplectic ODEs, it is beneficial to use 
symplectic integrators, due to the near-preservation of energy and excellent long-term beha- 
viour of the numerical solutions obtained by symplectic methods 1 7 , chapter VI] . Hamilton's 
principle states that the solution of a Lagrangian (in many cases also Hamiltonian) system 
moves along a path which extremizes the action integral S = L{q{t), q{t)) df among all 
paths q with fixed end points. 

One technique for deriving symplectic methods is based on the notion of discretizing 
Hamilton's principle, that is, replacing the action integral with a discrete action sum, and 
extremizing over all discrete paths or sequences of points qo,qi,...,q^ with fixed end points. 
These methods are known as variational methods or variational integrators. Variational meth- 
ods are guaranteed to be symplectic since the terms [qt-i, qic) of the discrete action sum can 
be interpreted as generating functions (of the first type) for the numerical flow map. Variational 
methods have been studied by numerous authors, we refer to the review article by Marsden 
and West 1 1 1 1 or the more recent encyclopedia article by Leok 1 10 1 and the references therein 
for more information about variational methods. 

Standard Lie group methods, like RKMK methods or CG methods, give numerical solutions 
that evolve on the same manifolds as the exact solutions. The question of the existence of 
symplectic methods of formats similar to the ones considered by Crouch and Grossman or by 
Munthe-Kaas has been a topic of interest for several years. 

On IR" there is a unified way to extend Runge-Kutta methods on U n to symplectic methods 
on T*U n |7, Section VI. 6.3], i.e. symplectic partitioned RK methods. Our goal with this article 
is to construct and study symplectic methods of arbitrarily high order that are extended from 
Lie group methods, i.e. high-order symplectic Lie group integrators. We have focused on the 
case where M = G and the action is simply multiplication in the Lie group. 

The idea of constructing variational methods from Lie group methods has previously been 
considered by several authors. Bou-Rabee and Marsden proposed in 2009 to base variational 
methods on RKMK methods |1|, and present a class of methods of first and second order. 
Methods of a similar type are also considered in the survey article by Celledoni, Marthinsen and 
Owren [3] . In the present article, this idea is pursued further to obtain methods of arbitrarily 
high order. 

It is known to the authors that a different, but related approach to variational Lie group 
methods has been studied by Leok and collaborators. The idea appears already in Leok's 
doctoral thesis |9, Section 5.3], and also in other articles by Leok and co-authors. A more 
detailed study of this approach will be in a forthcoming article by Hall and Leok. 
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1.2 ODEs on Lie groups 



Let G be a finite -dimensional Lie group and q its associated Lie algebra. We denote right- 
multiplication with g eG as Rg and left-multiplication with g as Lg. We use dot notation to 
denote translation in the tangent bundle, i.e. 

g ■ v = TLg v, v-g = TRg v, ve TG, 

and in the cotangent bundle 

g ■ p = T*L g -i p, p-g= T*R g -i p, pe T*G. 
We also need the notation Adg := TLg o TR g -i . All autonomous ODEs on G can be written as 

g = f(g)-g> g(0) = go, (1) 

where g is a curve in G and the map /: G — ► q is determined uniquely by the vector field. We 
can solve this kind of equation numerically using Lie group methods |8|. Here we have chosen 
the right- trivialized form of this equation. We could also have used the left-trivialized form 
g = g • f[g) which would have resulted in only minor changes to the formulae presented later 
in the article. 

Since we are interested in solving Hamiltonian ODEs using Lie group methods, we need a 
group structure on the cotangent bundle of G, as well as the map / that corresponds to this 
type of ODEs. 

1.3 Hamilton-Pontryagin mechanics 

Lagrangian mechanics on G is formulated in terms of a Lagrangian L: TG -» R. Hamilton's 
principle states that the dynamics is given by the curve q: U —- G that extremizes the action 
integral 

S H = / Uq,q)dt, 
Jo 

where the endpoints g(0) and q[T) are kept fixed. In fil Theorem 3.4] it was shown that this is 
equivalent to the Hamilton-Pontryagin (HP) principle, which states that the dynamics is given 
by extremizing 

Shp = I {Uq, v) + (p,q- v)) dt, 
Jo 

where v e T q G, p e T*G are varied arbitrarily, and the endpoints of q are kept fixed. This action 
integral leads to dynamics formulated on T*G. To simplify calculations, it is convenient to 
right-trivialize T*G to G x q* via the map {q, p q ) <-* {q, p q ■ q~ l ). Letting £(q,0 := L{q,£ ■ q), 
£, e g, it is easy to show that the HP principle is equivalent to the right-trivialized HP principle, 
which has action integral 

S= [ T [£{q,0 + (p,q-q~ l -0)dt, 
Jo 
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where ^: U —> q and p: U -* q* are varied arbitrarily, and the endpoints of q are kept fixed. Taking 
the variation of S, we arrive at the right-trivialized HP equations 

q = $-q, 

/i=-ad|/i + (Di^,f))-<y _1 , (2) 
p = D 2 £(q,S), 

where D k £ denotes the partial derivative of £ with respect to the fcth variable, i.e. a one-form. 
This is the ODE onGxg* that we need to solve. 

1.4 Variational integrators 

Variational integrators are constructed by discretizing an action integral and then performing 
extremization with fixed endpoints. This procedure turns the action integral into an action 
sum. The discrete Lagrangian L^: G x G — U is an approximation of the action integral over a 
small time step h, 

nkh 

Lh(qk-i><lk) ~ / L{q,q)dt, 
J(k-i)h 

where q: R — G extremizes the action integral with q(Q) and q(T) fixed. Letting N= Tlh, the 
action sum becomes 

N 

Sh = Y, L h(qk-i>q0- 

k=l 

Extremizing Sh while keeping qo and q^ fixed gives us the discrete Euler-Lagrange equations 

Y> x L h {q k , q k+l ) + D 2 L h (q k - l , q k ) = 0, 1 < k < N. 

The discrete Legendre transforms define the discrete conjugate momenta 

Pk '■= Pk ■ Ik '■= -ViL h {q k , q k+ i), 
Pk+i '■= Pk+i • Qk+i '■= V 2 L h {q k , q k+ i). 

By demanding that these two definitions are consistent, we automatically satisfy the discrete 
Euler-Lagrange equations. If we can solve the first equation for q k+ \, we can use the second 
one to calculate jUjfc+i, giving us the variational integrator iq k ,p k ) >-» {q k+ i,[i k +i). 

2 From a Lie group method to a variational integrator on the 
cotangent bundle 

2.1 Group structure and Hamiltonian ODEs on Gxg* 

Now that we have the ODE that we need to solve, we can choose a group structure on G x g* 
so that the map / takes a suitable form. We will start by assigning the following group product 
to T*G: 

{g, p g ){h, p h ) := {gh, pg-h + g-p h ). (3) 
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We note that 



(g, p g ) (h, p h ) = (gh,pg-h + g- pi 

= (gh,[p g -g- 1 +Ad* g _dPh-h- l ))-gh 

Thus, letting p = p g - g -1 and y=p\ l -h~ y , the right- trivialized version of J3) is the product on 
G x g* defined by 

(g,/i)(fe,v) := [gh,p+ Ad*. 1 v). 

It can be shown that the Lie algebra belonging to the Lie group Gxg* is gxg* equipped with 
the Lie bracket [[£,p), (ry, v)] = (ad^ry.ad* /U-ad^ v), where ad x is the derivative of Ad exp (x) with 
respect to x at the origin. We will also need an expression for TR Z ( for z = (q, p) e G x g* and 
(={T],v)eQxg*: 

TR Z (= ^[exp{eT]),ev){q,p) 



e=0 



= ^(exp(c77)^ev+Ad* xp( _ £)7) p) 
= [rj -q,v- ad* p), 



e=0 



We would now like to use this to write the right-trivialized HP equations J2) in the form of 
{]}, z = f{z) ■ z. If the map /: G x q* —■ q x q* satisfies 

f{q,D 2 £{q,0) = [Hi£(.q,Qi) ■ q~ l ) (4) 

for all {q, £,) e G x g, we see that z = f{z) ■ z, which is exactly what we need. 

In many cases, the map [q,0 >-* (q,~D2£{q,0) is a diffeomorphism of manifolds. If this 
holds, we say that the Lagrangian £ is regular. If £ is regular, the Lagrangian problem has an 
equivalent formulation as a Hamiltonian ODE on T*G, where 

JC(q,n 2 £{q,0) = (D 2 £{q,0,0-£iq,0 

and 

f{q, p) = \p 2 Je{q, p),~ [DiJ?(q, p))^' 1 )- 

For Hamiltonians which arise in this manner, the map {q,p) >-» (q,D 2 <ffl[q,p)) is also a dif- 
feomorphism of manifolds. In fact the map is the inverse of the one above. Hamiltonians for 
which this hold are also called regular. 

The group structure in (J3j was used by Engo in |5| to construct partitioned Runge-Kutta- 
Munthe-Kaas methods on T*G, without any special regard to symplecticity. 



2.2 General format for our integrators 

It is natural to consider discrete Lagrangians based on approximation of the action integral by 
quadrature. Consider the discrete Lagrangian 

L h (q ,q l ) = L h (Q l ,...,Q s ,S l ,...,S s ) = h£ i b i £(Q i ,Si), 

i=l 
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where Z?, are quadrature weights, and the auxiliary variables Qi,...,Q s ,£i,---,%s are chosen to 
extremize under the constraints 



Y(Q l ,...,Q s ^ 1 ,...^ s ,qo)-\og[q l q Q 1 ) = 0, 
X i (Qi,...,Q s ,Si,...,S s ,q Q )-log(Q i q Q ~ l ) = > i = l,...,s. 



(5) 



The functions Y and Xi will typically arise from Lie group integrators, as we will see later on. 
Let A be a Lagrange multiplier corresponding to the constraint containing Y, and let A ; - be 
Lagrange multiplier corresponding to the equation containing X, for i = 1, . .. , s. To obtain a 
variational integrator, we must extremize 



4-(yl,F-log^ 1 ^- 1 ))-X:(AoX-log(Q^o 1 )), 

i=l 

while keeping qo and q\ fixed. Varying this with respect to A, A;, and Qi, we obtain the set 
of equations 

qi = exp(Y)g , 
Qi = exp{Xi)q , 



dLh 

dLh 
dQi 



dY\* 



i ay 



A, 



(6) 



J 



dQ 



i I 



Ay-((dexp x J)*A/)-Q f> 



for all i=l,...,s. 

To define the integrator based on the generating function Lh, we need to evaluate the partial 
derivatives of Lh with respect to qo and q\. In doing so, we consider Qi,...,Q s ,£i,...,£s as 
functions of qo and q\ defined implicitly by |5j and |6j. The partial derivatives of Lh are then 



dL h Q Qj dL h d{j 



dq j\dQj 9^0 3£/ dq 



dL 



az^ aQ ; - a4 a^ 



(7) 



a^i ^ iaQ ; - 8(5,1 9( v/' ^ 



The functions Qi,...,Q s ,f satisfy the constraints {5} for all qo, q\. By differentiating the 

constraints we see that the identities 



dY QQj ar a^ y 

[dQj dq d£j dq 



+ dexp_ l Y oTR qo i, 



dY „ 

dq t 
_faF 3Q; ar a^\ 

= E ^r°^ + ^r°T- L "dexpy l °TR-i, 
jVdQj dqi a^- a^J * 



ex,- ^ 



'ax 3Qy ax,- a^ 

3^o Y \dQj dq d$j dq , 
(dXi dQj dXi d$j 



+ dexp_^- o TR a -i - dexp^ 1 o TR -i o 



3Qi 



(8) 



= V — -o— ±L + — -o-^- -dexpy X ori? n -io--^ 
j[dQj dq, dtj d qi ) yx - Q - dq, 



dQi 



v < dq 



i = l,...,s, 
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all hold. 

We combine the discrete Legendre transforms 

with l[6) and J7J, and simplify using {8) to obtain the equations 

Mo = |(^) ^ + % •'?o 1 + (dexp: 1 y )*yl + X(dexp:^)*A j) 

jUi = (dexpy 1 )*^. 
Using the identity J|), we get 

f{Q i ,B 2 £(Q i ,t i )) = (z i ,lB l £[Q i ,{ i ))-Q7 l) j, 
and defining m, Mi e g* by 

^ = hbiD^iQutt) = hbifii ■ Qi, ^ = hbiD 2 £iQi,$0 = hbiMu 



for i = l,...,s, we get 



97 V 



i I 



•(^-(dexp^Aj, 



(dY\* 



if ^1 



7 V a^ f 



Combining everything above, the variational integrator is defined by the set of equations 



Mo 



hbim = 



hbiMi = 



dq j j\dqol 
i dY\* 



(dY\ 



(dXj\ 



* \ 



A 



• % 1 + ( dex P- y) * A + Z ( dex P-x,- ) * x i > 
i 

•QrMdexp-^A*, 



(9) 



Q ; - = exp(X;)^ , i = l,...,5, 

<7i = exp(F)^ , 
jUi = (dexp^ 1 )*^. 

Notice that we no longer involve the Lagrangian. We only need to evaluate the vector field 
through the map /. This opens up the possibility of applying the method to Hamiltonian 
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systems which do not have an equivalent Lagrangian formulation (or indeed to any ODE on 
T*G). 

It should be noted that since the integrator can be formulated as a variational integrator 
on G, the group structure chosen for T*G in (3) is not consequential. For any choice of group 
structure on T*G such that the canonical projection T*G —- G is a homomorphism of Lie 
groups, there is an equivalent formulation of the integrator in |9) . Note that / is defined via 
the group structure, and a change of group structure would lead to / being changed as well. 

3 First and second order variational integrators 

In the article by Celledoni et al. (3) , a special case of variational integrators of the form in- 
troduced in the previous section was considered. These integrators serve as an example of 
application of the formulae above. In these methods, let aij and bi be the coefficients of an 
5-stage Runge-Kutta method, satisfying fo, ^ for all i. Let the discrete Lagrangian be given by 

L h {q ,q l ) = hY j b i £{Qi,^), 

i 

and the constraints by {5} and 

Y=hY j b i $ i , X i = hY j a i j^ j , i = l,...,s. 

i j 

We can see that for i,j = l s, 



dY 
dq 


= 0, 


djq 
dqo 


= o, 


dY 
QQi 


= 0, 


dXj 
dQi 


= 0, 


dY 

Wi 


= hb u 


dXj 


= haji 



By inserting these into J9), we get the set of equations 

Mo = (dexply)*/l + ^(dexp:y *A ; -, 

hbitii - -(dexp^)*A;, 

hbiMi = hbiA + YjhajiAj, 
j 

{$i,n i ) = f{Q i ,M i ), 

Qi = exp{Xi)q , i = l,...,s, 

q-L = exp{Y)q , 
jUi = (dexpy 1 )*/!. 
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In these equations, A and Xj can be eliminated, the integrator can be written as 

biMi = bidexp*_ Y L + hj^bj Ad* emXj) nj) - h^bjajidexp*^ n } , 

j j 

j 

Qi = exp(Xi)q , 
tfi,ni) = f(Qi,Mi), i = l,...,s, 

Y=hY,bj{j, 

i 

qi = exp(F)g , 

Mi = Ad* xp( _ Y) [^o + h £ bj Ad*. 



(10) 



exp(Z,) n j 



Here we have used the identity dexp x odexpl!,. = Ad exp (x)- Equation (To) is equivalent to the 
method presented in (3j Section 5] . 

Methods of this form suffer from an order barrier. They can not obtain higher accuracy 
than second order. The proof, presented below, is closely related to a similar order barrier 
for commutator- free Lie algebra methods |2|. In the proof, a non-regular Hamiltonian is 
introduced. The Hamiltonian ODE does therefore not have a corresponding Lagrangian formu- 
lation, but as noted earlier, the methods derived in this article can be applied to it nevertheless. 

Proposition 3.1. The integrators of the format i fTo) can not achieve higher than second order 
on general Hamiltonian differential equations. 

Proof. Let a Hamiltonian on G x q* be given by J€{q,[£) - (/j,, v{q)), where v. G —- q is smooth, 
but otherwise arbitrary. The corresponding Hamiltonian vector field is 



f{q,H) = [v{q) t - 
and the differential equation is 



(dv_ 
dq 



q = v{q)-q, 



[dq 



If [q[t),iJ.(t)) is a solution to this ODE, then q(t) is independent of /Li, and furthermore, q[t) 
solves the ODE q = v{q) ■ q. Applying |To] to this particular Hamiltonian equation yields a 
numerical solution where q\ is independent of jUq and is defined by the equations 



£/ = v{Qi), 

Qi = exp\hJ^aijZj}q , 

j 

qi = exp\hj^biti}q . 



i = l,...,s, 
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We recognise these equations as a commutator-free Lie group method with one exponential, 
or equivalently, an RKMK method with cut-off parameter 0, applied to the ODE q = v[q) • q. 
As explained in 1 14|, commutator- free methods with one exponential cannot satisfy the third 
order conditions, and the solution is at most second order accurate. 

By applying the variational method (To) to a particular class of Hamiltonian problems, and 
considering only the Lie group part, or the projection, of the solution, we have an error of at 
most second order, thus the variational method is at most of second order as well. □ 

By repeating the argument with any variational integrator of the form described in J9) , we 
get a generalization. 

Proposition 3.2. A variational integrator on T*G based on a Lie group integrator can not 
achieve higher order than the Lie group integrator on general Hamiltonian differential equa- 
tions. 

The proof of Proposition|3. l|on the preceding page uses a Hamiltonian which is not regular 



and which does not have a corresponding Lagrangian formulation. Proposition 3.1 can also 
be shown to hold for regular Hamiltonians by considering a right- invariant Lagrangian, i.e. a 
Lie-Poisson equation, and expanding numerical solution pi and exact solution p[h) in Taylor 
series of h. This expansion is straight forward, but tedious, and is not included in this article. 
The proof of the more general Proposition |3.2| uses the non-regular Hamiltonian. It seems 



likely that Proposition 3.2 holds for the case of regular Hamiltonians as well, but we do not 
have a general proof of this. 

We should note that first and second order methods of the format described in fTol do exist. 



Specifically, a method of that format is first order if hi = 1 and second order if, in addition, 
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bittij = 1/2. The proof is a special case of Theorem 5.2 on page 

Example 3.3 (The midpoint method). Let us choose s = 1. The method is of second order 
if and only if we choose h\ - 1 and a u = 1/2. This method is also symmetric, since the 
substitution h —- -h, qo *-* q\, po p\, Y —- -Y in flO) yields the same method after some 



manipulation of the equations. This property is utilized in Section 6.1 to achieve high order by 
composition. 

4 Higher order integrators 

The methods discussed in the previous section were limited to at most second order. To obtain 
higher order integrators, we consider two approaches, based on two well-known classes of Lie 
group integrators. 

The first approach is based on the Runge-Kutta-Munthe-Kaas (RKMK) methods. This 
approach was already considered by Bou-Rabee and Marsden in 2009 1 1 1 . The work presented 
in the present article builds on the work by Bou-Rabee and Marsden, and examines in detail 
the case when the cut-off parameter r [q in |T] 12 ) in the RKMK method is larger than 0, and 



provides a complete order theory for variational methods based on RKMK methods. 

The second approach is based on Crouch-Grossman (CG) methods. This approach has to 
our knowledge not been explored before. We show that these methods can achieve arbitrarily 
high order, but the complete order theory of these methods remains unresolved. 
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4.1 Variational Runge-Kutta-Munthe-Kaas integrators 

A popular class of Lie group integrators is the class of RKMK integrators. For our purpose, these 
integrators can be written 

Qi = exp(Xi)q , 

ti = f(Qi), i = l,...,s, (11) 

Y=hY J bidexpJ r ) Xi Z i , 

i 

qi = exp(F)^ - 

where 

1 r Bjc 

dex P«,x = id - - ad x + £ "jy (ad*) k 

is the Taylor series approximation to dexp" 1 , and aij,bi are the coefficients of a Runge-Kutta 
method. If the RK method is of order p and r>p-2, the resulting Lie group integrator is of 
order p as well |7, Theorem IV8.4]. 

Variational methods based on RKMK methods were considered by Bou-Rabee and Marsden 
in fil, though the methods they present in detail are at most second order, since they only 
consider the case r = 0. The methods in this case are essentially the methods considered 
in Section|3j For r > 0, some complications arise for the variational integrator, since the X\ 
are not explicitly given by £i , . . . , f s . Our solution is to treat both x ; and £ , as unknowns and 
the equations for xi in l[TT} as restrictions. The Lagrange multipliers A,- corresponding to 
the equations for xi cannot be eliminated from the equations in a general manner, so the 
dimension of the non-linear equation to be solved at each step is larger than that of the 
corresponding symplectic method applied to T* U n , or the simpler integrator with r = 0. 

In the variational integrator we set the discrete Lagrangian to 

L h = hY 1 b ] £{Q ] ,$j), 

i 

and let the constraints be given by l[5j, that is 

qi = exp(F)^ - 
Qi = exp{Xi)q , 

and 

Y = hY,b J dexp- ( - r ) >Xj S j , 
i 

X i = hY f ciij dexp"} x . i = l,...,s, 

where X{ - log[Qiq^ l ). Note that on the solution set of the constraints, X; = jc,-. 
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Applying the variational equations from (9), the integrator is given by 
fJ-o = ((dexpl^)* - fc£>(dexp:y * oP^XtJi^A 

i 

+ l((de J q ) :y-ft£fl // (dap:yoP ( * r) (X / ,w)Aj, 



n; = - (dexp z |) A; + fefo; (dexp^- 1 ) o P* ) {Xi,^)A 
+ hY^ajildexpx])* °P* r ) CX/>fr)A;, 



(12) 



^fciM; = h(dexp (r 5 X! )*(^yl + ^a 7i A 7 -J, 

Q; = exp(X / )^ , 
fo,ni) = /(Q/,Mi), i = l,...,s, 

= exp(F)g > 
Hi = (dexpy 1 )*^, 

where P* r) (x, if) is a polynomial in ad* and ad^ of degree r, defined as the adjoint of the partial 
derivative of dexp^ x £ with respect to x. Specifically, 

P* 0) (x,$)=0, 

P* 2) = 1 ad* - - ad* ad* x + — ad* ad* , 
^^a^adl-E^Ead* (ad*)^- 1 . 

z k=2 K - i=0 

By applying Ad* xp(X ) to both sides of the second equation in (12) , we see that the first equation 
can be simplified. Using this and rearranging the rest of the equations, we arrive at the set of 
equations 

A = dexp! F (^ + fc£ bi Ad* xp(Xi) «,), 
i 

Xi = - hbi dexp* Xi ni + hP* {r) {Xi,^i)\biA + Y,aji^i], 

i 

Mi = yM e ^(ixX( biA + X< a j iA i)' 

Xi^hY^aijdexp^x.Sj, (13) 
i 

tfi.rii) = f{exp{Xi)q ,Mi), i = l,..,,s, 

Y = hY J b i dexpJ r ) Xi ti, 

i 

qi = exp{Y)q , 

Ml = Ad* xp( _ y) (MO + b i Ad exp(Jf,) n i)> 
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which define a symplectic integrator on T*G. The identity dexp x o dexplj. = Ad exp ( X ) was used 
to obtain the last line of the equations above. We call the integrators defined by |T3| variational 
Runge-Kutta-Munthe-Kaas methods or VRKMK methods for short. One can easily check that if 
the Lie group is abelian, a VRKMK method simplifies to a symplectic partitioned Runge-Kutta 
method. 

In implementations of these methods, it is required that exp: g -* G, dexp* : gxg* — ► g* and 
Ad* : G x g* —- g* are calculated to machine precision to obtain symplecticity. The equations 
defining the integrator can be solved as a set of non-linear equations in the unknowns X, , 
and Ai, i = 1 5, as the other quantities in fl3) are given explicitly in terms of the afore- 
mentioned variables as well as qo and hq. If G is an n-dimensional Lie group, this is a total of 
3ns scalar unknowns. Other choices of independent and dependent unknowns are possible. 

The number of unknowns can be reduced by using a (1, r)-Pade approximation of dexp , 
that is dexp~j , = l/(l + ^adf h — ), instead of the Taylor series approximation above. The 
details of this alternative approach is not included in this article, but should be explored if 
efficient implementation is an objective. Alternatively, one could obtain explicit expressions 
for Xi and therefore be able to eliminate the X\ in the integrator by starting with coefficients of 
an explicit RK method. The variational method based on an explicit RK method would still be 
implicit, however, and due to the order conditions presented later in this article, the increased 
number of stages required for a particular order would offset the reduction in number of 
unknowns per stage obtained by using an explicit method as the underlying method, so it is 
unclear if this simplification is useful in practice. 

4.2 Variational Crouch-Grossman integrators 

The methods of Crouch and Grossman form another important class of Lie group meth- 
ods. Crouch and Grossman formulated their integrators in terms of rigid frames, i.e. finite 
collections of vector fields on a manifold. On a Lie group, a suitable rigid frame is a basis 
for the right-invariant vector fields on G corresponding to a basis of g. In this setting, the 
Crouch-Grossman methods can be defined as follows. Let bi, a\j be the coefficients of an 
5-stage RK method. The Crouch-Grossman method with the same coefficients is defined by 
the equations |7, Section IV 8.1] 



Qi = exp(ha is Ss) 

ti = f<Qi), 

qi = expihbsts)- 



■exp(han(i)q , 



exp{hb^i)q . 



The order of a CG method is determined by the order conditions developed in 15 . 



Using the general format J9) , we set the discrete Lagrangian to 



Lhiqo, qi) = h £ &^(Qi,£i), 



(14) 



with bj ^ and constraints given by JH), that is 



qi = exp(F)^ , 
Qi = exp[Xi)q , 
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and 

Y = log[exp{hb s Z s ) ■ ■ ■ exp{hbiti)), 
X i = log(exp{ha is Z s )---exp{hanZ l )). 

Inserting this into the equations defining a variational integrator J9), we obtain 

,uo = (dexpI 1 y )*yl + ^(dexp:^.)*A ; -, 

j 

hbirii = -(dexp^)* A,-, 

hbtMi = hbi dexp*^. °Ad* exp(hb ^ sh . exp{hbi+i(i+i) o(dexpy l )*A 

+ aji dexp* a . jfj °Ad* expihajs(s) ... exp[haj °(dexp^)* Ay, 

Gi,ni) = f(Qi,Mi), 
Qi = exp{Xi)q , 
qi = exp{Y)q , 
Hi = (dexpy 1 )* A 

Eliminating A and Ay, and rearranging, we get the integrator 

q x = q s , q j = exp(Jibjtj)q j ~ l , q° = q , 
Qi = Qis, Qij = exp{ha i jZ j )Q i j- 1 , Qi = q , 
tfi,ni) = f<.Qt,Mi), 

Ho = Ad* o Ho, Hi = Ad * qi Ml. m = Ad* Qi tit, 

s 

Hi = fia + hY^bjfij, 

7=1 



Mi 



y_2 U l 



(15) 



for all i = 1, ... , 5. We call the integrators defined by \T5f variational Crouch-Grossman methods 
or VCG methods for short. The last equation in {15) can also be written as 



In the case that the Lie group is abelian, the integrator simplifies to the same symplectic, 
partitioned RK method as in the abelian case for the VRKMK integrator. 



5 Order analysis 

In analysing the order of variational methods we use variational error analysis as described 
by Marsden and West 111, Section 2.3] . We recite two definitions, and a theorem from the 
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reference above which are useful in the following sections. A discrete Lagrangian is said to 
be of order p if 



L h [q{0),q{h)) = [ L{q{t),qit))dt + 0(hP +1 ), 
Jo 



for all solutions q{t) of the Lagrangian problem. Two discrete Lagrangians L l h and L 2 h are 
equivalent if their difference, L l h [qo, q\) - L 2 h [qo, q{], is a function of h only, i.e. independent of 
qo and q\ . If two discrete Lagrangians are equivalent, the symplectic mappings based on them 
are identical. 

Theorem 5.1. Given a regular Lagrangian L and corresponding Hamiltonian H, the following 
are equivalent for a discrete Lagrangian L^: 

• The symplectic integrator defined by L^ is of order p. 

• Lh is equivalent to a discrete Lagrangian of order p. 

Both classes of methods presented in this article depend on Butcher coefficients a,-; and £>,-. 
Furthermore, when applied to an abelian Lie group (for instance U"), both classes become 
symplectic partitioned RK methods where the position is integrated with the RK method 
with coefficients aij and hi, while the momentum is integrated with a method with coeffi- 
cients dij = bj - bjaji/bi and hi = hi. These methods can be called variational Runge-Kutta 
methods, and the order conditions for these have been explored in detail by Murua 1 13 1 . Since 
an abelian Lie group is a special case of a Lie group, the order of the variational RK method is an 
upper bound for the order of the variational Lie group method with the same coefficients. The 



order of the underlying Lie group method is also an upper bound according to Proposition 3.2 
on pagefTQ] 

5.1 Order of VRKMK integrators 



The methods described in Section 4.1 are fully described by the Butcher coefficients aij and bj, 
and the cut-off parameter r. The cut-off parameter r limits the order of the RKMK method fTT) 
on G. The order of the RKMK method is the minimum of the order of the RK method based on 
the same coefficients and r + 2. |7, Section IV.8.2] 

As explained above, the order of the VRKMK method is bounded from above by the order of 
the variational RK method based on the same Butcher coefficients, and by the order of the 
RKMK method. Since the order conditions for RK methods for a particular order form a subset 
of the order conditions for the variational RK method, this means that we can a priori say that 
the order of the VRKMK method is bounded from above by the order of the variational RK 
method and r + 2. The following theorem states that the order of the VRKMK method is in fact 
the minimum of these two bounds. 

Theorem 5.2. If the variational Runge-Kutta method based on the coefficients aij and hi is 
of order at least p, and the cut-off parameter r satisfies r>p-2, then the variational Runge- 
Kutta-Munthe-Kaas method with the same coefficients is at least of order p. 
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Proof. The proof consists of two steps. In the first part, we show that the discrete action sum 
arising from the "full" RKMK method without cut-off, that is, the method where dexpj" * x is 
replaced by dexp" 1 , approximates the action integral to order p, ensuring that the correspond- 
ing variational method is of order p. This is done by considering a Lagrangian problem on 
Tq, whose action integral is estimated by the underlying RK method. In the second part, we 
show that the difference between the VRKMK method with cut-off parameter r and the "full" 
VRKMK method goes to zero as 6{h r+2 ). To distinguish between the two RKMK methods, we 
will denote the "full" RKMK method RKMK(oo), and the RKMK method with cut-off parameter 
RKMK(r). The variational integrators based on the two methods are denoted VRKMK(oo) and 
VRKMK(r), respectively. 

Let q: [0, a] — G be a solution to the Euler-Lagrange equation with q{0) = qo, and assume 
that a > is sufficiently small such that a{t) = log(<7(t)«7 ( ^" 1 ) is uniquely defined for all t e [0, a]. 
The exact discrete Lagrangian is given by 

L E h {q ,q{h)) = [ k £{q{t),at))dt, (17) 
Jo 

where ^(f) = q{t) ■ q{t)~ l . If we define h Tg — R as 

£[a,&) = £{exp{a)q ,dexp a &), (18) 

we can rewrite (17) as L^(qo, q{h)) = Z^(0,(7(h)) where 

r h 

L E h {0,a(h))= i(aW,a{t))dt. 
Jo 

This is an exact discrete Lagrangian on the vector space g, which we approximate by the action 
sum arising from the underlying RK method. 

l^{OMh)) = hY i biiiy i ,TiO, (19) 

where y; = hY.jO,ij'qj> i = 1 5 and the sum is extremized under the constraint a{h) = 

hHi biT]i. Under the assumptions of the theorem, the order of the variational Runge-Kutta 
method is at least p, so 

If {0Mh)) = Ll[o,a{h)) + ©{HP +l ). 

Inserting (18) into fl9) gives 

Lf[q ,q(h)) = hY i b i £(exp(.y i )q ,dexp yt i]i), (20) 

i 

where the sum is extremized under the constraint q{h) = exp[h JJ,- birji)qQ. 
Now, the discrete action sum arising from RKMK(oo) is 

Lf MUoo \q Q ,q{h)) = hY,bi£{exp(Xi)q ,$i), (21) 
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which is extremized under the constraints 

X i = hY J a-ij dexp^ 1 £/, i = l,...,s, 



q{h) = exp^hj^ bi dexp^ 1 £jj q . 



We see that under the identifications 

yi = Xi, r]i = dexp^* 

the objective functions (20) and |2~fl are identical and are extremized under the same con- 
straints. Thus their extremal values are identical and we have proved 



r RKMK(oo) 



{qo,qm) = L™[q ,q{h)) 



= Ll K {Q,a(h)) 

= L E h [0,a(h))+©{hP +l ) 

= L E h {q ,q{h))+6{ht> +1 ), 

concluding the first part of the proof. 

For the second part of the proof, we consider the integrator in (13} and the variational 
integrator based on RKMK(oo) with the same initial data [qo, /io). Let £/, ni,Xi, Mi, Xi,A and Y 
be as in fT3), and etc. be the corresponding quantities in VRKMK(oo). 

We define 8^ = - cf (oo) and so on, and consider the difference between VRKMK(oo) and 
VRKMK(r). Since q\ = exp{Y)qo and ^ = (dexpy 1 )*/!, the leading order of the difference 
between the two integrators is given by the leading orders of 8Y and 8 A. It is clear from the 
equations in fl3) defining the integrator that as h — 0, X it Xf^fX^X^, Y and F (oo) all go to 
zero as 6{h). Furthermore, 5^,5n;,5M,- and 8A must also go to zero as h —- 0. By using the 
expressions in fl3} and the series expansions of dexp* and Ad* xp(i) , we find that 



571 
SMi 
8Xi 



ad^jUo + h^bi[adg x . n,[ + 5n,) + higher order terms, 



. I ad * A _ fad* ) r+1 /i + 8A + Y ^8Xj + h.o.t., 



SXt 



---hbi 

+ hbii- 
+ h.o.t., 



hJ^Uij 
3 



aA* sx . rii + Srii 



(Xi,ti)A + 



-ader — adt adl Y H — ad^ y ad, 
U df( 6 (i dXi 12 6X1 (i , 



A + - ad? 5/1 

2 f < 



(r + D! 



(adi/ + %-fadJ x f; + «fj 



+ h.o.t., 



(fffi,ffn/) = T iQiM) f{8Xi ■ Q t ,SM0 + h.o.U 



SY = hY J b i 



B r +i 



(r + 1)! 



(r+U 



f/--adJ Jti ^ + fffi| + h.o.t, 
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where Q* r+1) ix, = P* r+ i) ( x > ~ P* r ) (x, £)• In the equations above, "higher order terms (h.o.t.)" 
denote terms that are dominated by at least one of the preceding terms on the right hand side. 
From the definition of P* r) [x,0 and X t = 0{h), we see that Q* r+1) = @{h T ). 

We know that 5^- and 8 ni go to zero at least as 0{h l ). We proceed by induction and assume 
that 8£i = 0(h k ) and Snt = 0(h k ) for all i. From the equations above, we can then successively 
see that 



and finally, 



SXi 


= 0{h k+1 ) 


>+0(h r+2 


8Y 


= 0[h k+1 ) 


i + 0(h r+2 


5A 


= 0[h k+1 ) 


i + 0(h r+2 


8Xi 


= 0[h k+1 ) 


\ + 0{h r+l 


SMi 


= 0[h k+1 ) 


\ + 0{h r+l 


SSi 


= 0[h k+1 ) 


\ + 0{h r+l 


8 m 


= 0[h k+1 ) 


\ + 0{h r+l 



By an inductive argument, we can conclude that 8$i = 0{h r+l ) and 8ni = 0{h r+l ). Therefore 
SY = 0{h r+2 ) and 8A = 0{h r+2 ), concluding the second part of the proof. □ 

An immediate consequence of the proof is that there exist methods in this class of arbitrarily 
high order. For instance, the Gauss methods \7, Section II. 1.3] form a class of Runge-Kutta 
methods which achieve arbitrarily high order. Since these methods themselves are symplectic, 
d(j = bj-bjUjilbi = dij, and the variational method based on a Gauss method is a partitioned 
Runge-Kutta method with the same coefficients for both position and momentum. The vari- 
ational method is equivalent to the Gauss method itself applied to the Hamiltonian ODE, and 
has therefore the same order as the Gauss method itself. When r is large enough, the VRKMK 
method based on the coefficients of a Gauss method achieves the same order as the Gauss 
method. 

5.2 Order of VCG integrators 

In this section, we will prove that there exist VCG integrators of any order. To show this, we will 
need the following lemma. 

Lemma 5.3 (Composition of VCG integrators). Let {A m , b m ) and {A {2) ,U 2) ) be the Butcher 
tableaux of Runge-Kutta methods with and s (2) stages, andy a real number. The composition 
method formed by first applying the VCG method based on , with step length yh and 
then the VCG method based on [A {2 \ £> (2) ) with step length (1 - j)h, is a VCG method with 
Butcher tableau 
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(l-y)A (2) 



jb\ 



(i) 



(2) 



ii-y)b 



(2) 

.(2) 



Proof. Consider the discrete Lagrangians corresponding to the two VCG integrators that are to 
be composed, 



with constraints 



and 



..in 



i=l 



Q f w = eip(fcaJW n ^)---eip(fca«fS lJ )*. 
^ = exp{hb%S%) ■ ■ ■ exp{hb^^)q , 



with constraints 



i=\ 



(2) 



(2)r(2)^i _f }, ^,(2^(2)1 

s C9S 4 o 



^ = exp^f^-.-exp^^o, 
as well as the composition discrete Lagrangian 

Lfiqo, qi) = L^ h {q , q) + L^_ y)h {q, q x ), 



where q is chosen so that i' c) is extremized. It was proved in 



integrator corresponding to is the composition method that results from composing 
the integrator corresponding to L ^ with the integrator corresponding to L~ Denote by 
(A (c) , fo (c) ) the Butcher tableau with s (c) = s' 1 ' + s (2) stages given above. Then 







Theorem 2.5.11 that the 



Lf{q ,qi) = h 



,.(2) 



i=l 



t (c) 



2 = 1 
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with constraints 



Qf = exp{ha^%)-exp(haf^f)q , 
q l = exp(hb%S%)---exp(hb?^)q . 



□ 



Proposition 5.4. There exist methods of any order among the VCG integrators. 

Proof. From 1 7 , Section II.4] , we know that if we compose a one-step method with itself using 
different step sizes, we can obtain arbitrarily high order, provided we choose the number 
of steps and the step sizes appropriately. Thus, we obtain VCG methods of any order by 
composition. □ 

Proposition 5.5. For VCG integrators, the order conditions for first and second order are the 
same as for the underlying Runge-Kutta method, i.e. 



I> = 1, 



and 



1 



Ci = 



where Ci = atj. 

Proof. We use variational order analysis, as presented in [TT] Section 2.3] . Let the exact discrete 
Lagrangian be denoted 

nh 

L\{q a ,q{h))= \ L(q(t),q{t))dt, where q = q(0). 
Jo 

The exact discrete Lagrangian can be expanded in powers of h: 



L B h [q ,q{h)) = Y, 



k=0 



k\ 



' d k 



L\{q ,q{h)) 



h=0 



From the right-trivialised HP equations j2) and (4), it is straight- forward to show that 

^-L{q{t),q{t)) = ^-{^,0 = </U> + <M> = </ 2 (z),/i(z)> + //z, ^-fi(z)\ 
at at \ at / 

where z = {q,y) and /(z) = (/i(z),/ 2 (z)). Thus, letting (£o. n ) = /(z ) = /(<7o,Mo)» and using 
e{q,0 = Uq,q),weget 

L\{q , q{h)) = h£(q ,S ) + y <n ,£o> + j^tfo ■ a o) + |^("o - ad* o ^)^J +0{h?). 

Similarly, we can expand the discrete Lagrangian in powers of h by using (14) together with 
Qi\h=o = (Jo and £,U = o = to- 

oo fok i jfc 

L h [q ,q[h)) = £ 



fc=0 *! Idfc* 
k I d k 



L h {q ,q{h)) 

h-- 

OO I.K / S \ 

k=0 kl \dh k i=l h=0 ) 
= h(Y J b i )£{q ,^) + ^-l2Y J b i ^-£(Q i ,^ 



+ ©ih 6 ). 



h=0J 
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By comparing equal powers of the two expansions, we see that the first order condition is 
£j b l : = 1, as in RK methods. The second term needs more work. We apply (4) together with 
m \h=o = n and Mi\ h=0 = fx and get 



h=0 



dQi \ I d£ 



d/i 
dQi 
dfr 



d/z 



ft=0 

d^ 
dfc 



h=0 



We calculate the derivatives of Q; and f j with respect to h using (15) : 



d^ 
dh 



dQi_ 
dh 

h=0 



h=0 



dq dh 



3/i dMi 

+ - — o 



h=0 dfi dh 



h=0) 



We also need the derivative of . In this expression, we apply (TBJ and simplify using the first 
order condition: 

dMi 



dh 



h=0 



(l n "^ad* o fi . 



Putting these equations together, we obtain 



2 1> — ^(Q; , W = 2 £ hi a ( n , < > 
/ dn fe=0 Y 



+ ( no, 2Y biCi— -(£ ■ q ) + ^ 2f 1 - Y biCi)n - adt /z 



9/i 



Thus, to get second order, we need the second order RK condition = 1/2. 



□ 



The computation for third order is similar, but much more complicated. In general, the 
order theory for VCG integrators is not complete and needs further study. 



6 Numerical tests 

The methods were tested on the test problem "dipole on a stick" (see Figure[TJ, a pendulum 
consisting of a long straight rod of length 1, with one end fixed, but freely rotating, at the origin, 
and a shorter rod of length 2a attached perpendicularly to the long rod at the other end. At 
each of the endpoints of the shorter rod there are charged particles with masses ml 2 and 
electric charges ±q. The rods are assumed to be massless. The pendulum is affected by gravity 
in the negative e3 -direction and the electric field generated by a charged particle at position 
z = (0, 0, -3/2) T of charge /3. The physical constants for specific gravity and electric force are 
set equal to 1. 

If we let y+{t),y-{t) denote the positions of the positive and negative charge, respectively, 
the position of the pendulum can be described uniquely by the matrix g{t) e SO(3) such that 
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Figure 1: Dipole on a stick 

y±(t) = g(f)y+, where y® = (0, ±a, -1) T is the position of the two charged particles in reference 
or body coordinates. Using the standard identification of so (3) with IR 3 and of so (3)* with IR 3 
via the standard inner product, the state of the system [g, /i), can be represented with g e SO (3) 
as a 3 x 3 real matrix, and ^ e so (3)* as a vector in IR 3 . 
The right-trivialized Hamiltonian of this system is 

^(g>M) = \v T gI~ l g T V+ mejge 3 + qP{\\gy° + - zf 1 - \\gy°. - z\\~ l ), 

where J = mdiag(l + a 2 , 1, a 2 ) is the inertia tensor of the pendulum. 

6.1 Order tests 

The VRKMK methods that were tested are based on the 1-, 2- and 3-stage Gauss methods, and 
Kutta's third order method. These methods are defined by the Butcher tableaux and cut-off 
parameters in Table[TJ These methods can be shown to satisfy the extra order conditions for 
variational integrators to their respective orders. 

The order of the VCG methods were also tested. To obtain higher order, symmetric com- 
position of the midpoint method (which is symmetric) as described in |7. Section V3.2] was 
used. The Butcher tableaux of the resulting methods are the same as those of the fourth and 
sixth order diagonally implicit Runge-Kutta methods (DIRK) shown in Table|2j The paramet- 
ers yi , . . . , 74 were derived by Yoshida 1 16) . 

The methods were implemented in Matlab, using a modified version of the DiffMan 
package 1 6 1 for defining Lie algebra and Lie group classes and functions on these spaces. The 
sets of non-linear equations (13} and fl5) were solved by fix-point iteration. 
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(b) Kutta's third order method 
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r = 2 

(c) Fourth order Gauss method 



r = 4 

(d) Sixth order Gauss method 



Table 1: Butcher- tableaux of the RKMK methods tested 



In these tests we have used the data 



1 = 


p= 


1, 


0.1, 






1 














-1 





1 






/i(0) = g(0)/g(0) T e 2 . 

The initial data /i(0) is chosen so that the first component of /(g(0),^(0)) is &2- 

The errors in ;u(0.5) and g(0.5) with respect to a reference solution are shown in Figure[2j 
The errors plotted are || — Mref II 2 + Ilg _ grefll2> where the first IHI2 is the Euclidean vector norm, 
and the second is the subordinate matrix norm. The reference solution was calculated using 
the sixth order VRKMK method with step size h = 10~ 3 . The dashed lines are reference lines 
for the appropriate orders and are the same lines in the two plots. As is evident from the plots, 
errors from fixed point iteration dominates the errors for the 6th order methods when h is 
small, and the methods appear to obtain their theoretical order. 

Analytically, the second order VRKMK and VCG methods are actually identical. The im- 
plementations of the two methods are different, as the non-linear equations are set up in 
slightly different manners. The result of this is that the numerical solutions differ slightly. For 
this numerical test, the VRKMK methods are superior to the VCG methods, since the error 
constants are lower. 
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(a) Second order midpoint method 
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(b) Fourth order DIRK method based on triple jump 
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74 
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7i 



71 = 0.78451361047755726381949763, y 2 = 0.23557321335935813368479318, 
73 = -1.17767998417887100694641568, y 4 = 1.31518632068391121888424973 

(c) Sixth order DIRK method 

Table 2: Butcher-tableaux of the VCG methods tested 
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Step size Step size 

(a) Absolute error of VRKMK methods (b) Absolute error of VCG methods 



Figure 2: Order plot. Dashed lines are reference lines for the appropriate orders 
6.2 Long time behaviour 

The long term behaviour of the methods was also investigated. In Figure|3] the energy error of 
the numerical solution is plotted over the time span (0, 1000). We have used step size h = 0.01 
(10 5 integration steps). Only the second and fourth order methods were tested on this time 
span. As can be seen from the plots, the energy error is small, approximately 10" 3 for both 
second order methods, and approximately 10~ 7 for the fourth order VRKMK method and about 
10~ 5 for the fourth order VCG method. The energy error is about 100 times lower for the fourth 
order VRKMK method than for the fourth order VCG method, confirming the superiority of 
the VRKMK methods for this test problem. 

7 Conclusion and future work 

In this article, a set of equations defining symplectic integrators for ODEs on T*G were presen- 
ted, as well as two classes of integrators using these equations. The integrators obtained are 
formulated intrinsically on T*G, and any drift away from the manifold in numerical solu- 
tions is due to round-off errors. The integrators were developed as variational methods for 
Lagrangian problems, and are therefore symplectic when applied to Hamiltonian differential 
equations. Both classes that were studied, were shown to contain methods of arbitrarily high 
order, although the computational cost per time step increases with the order. Effective imple- 
mentation of the methods has not been a major goal in this article, we have instead focused 
on the properties of these methods. 

The two classes of symplectic methods are based on, respectively, the Runge-Kutta-Munthe- 
Kaas methods, and the Crouch-Grossman methods. The methods have a partitioned structure 
where the position on the Lie group is integrated by the Lie group method while the mo- 
mentum is integrated by formulae which involve various functions on g* . We can therefore 
say that these methods are partitioned Lie group methods and are Lie group methods in a 
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wide understanding of that term. To the knowledge of the authors, this is the first time that 
symplectic Lie group methods have been presented and studied in the level of detail done in 
this article. 

The reformulation of RKMK methods with a Pade approximation of dexp -1 is briefly dis- 
cussed in Section 4.1 This reformulation makes it possible to eliminate the Lagrange mul- 
tipliers Xt, which is beneficial for computational efficiency. Implementation and study of 
this approach would make the variational RKMK methods more competitive in terms of 
computational cost. 

Another class of Lie group methods is formed by the commutator-free Lie group methods. 
The approach described in this article can easily be used to formulate symplectic Lie group 
methods based on commutator-free methods. Formulation, implementation and study of 
variational commutator-free methods are aspects that can be pursued in the future. 
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